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Abstract 

Ordinary differential equations obtained as limits of Markov processes appear in 
many settings. They may arise by scaling large systems, or by averaging rapidly 
fluctuating systems, or in systems involving multiple time-scales, by a combination of 
the two. Motivated by models with multiple time-scales arising in systems biology, we 
present a general approach to proving a central limit theorem capturing the fluctuations 
of the original model around the deterministic limit. The central limit theorem provides 
a method for deriving an appropriate diffusion (Langevin) approximation. 
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1 Introduction 



There are two classical kinds of Gaussian limit theorems associated with continuous time 
Markov chains as well as more general Markov processes. The first of these considers a 
sequence {X N } of Markov chains that converges to a deterministic function X and gives 
a limit for the rescaled deviations U N = r N (X N — X). (See, for example, Kurtz (1971 



1977/78); van Kampen (1961).) The second considers an ergodic Markov process Y with 
stationary distribution ir and gives a limit for 



Z N {t) 



Nt 



h(Y(s))ds 



N I h(Y(Ns))ds } 
o 



for h satisfying J hdir = 0. (See, for example, Bhattacharya (1982) for a general result of 
this type.) 

There are many proofs for theorems like these. In particular, results of both types can 



be proved using the martingale central limit theorem (Theorem A.l). For example, in the 
first case, there is typically a sequence of functions F N such that 



M N (t) = X N (t) - X N (0) - f F N (X N ( 

Jo 



s))ds 



is a martingale, F N — > F, X = F(X), and F 



N 'x N ) 



F(x) « VF(x)(x 



N 



X 



for x 



N 



converging to x. If the martingale central limit theorem gives r^M N =>• M and U (0) 
i7(0), then (ignoring technicalities) U N should converge to the solution of 



U(t) = {7(0) + M(t) 



VF(X(s))U{s)ds. 



(1.1) 



In the second case, the assumption that J hdir = suggests that there should be a 
solution of the Poisson equation Af = —h, where A is the generator for Y, and then 



Z N (t) 



(f(Y(Nt)) - f(Y{0)) 



Nt 



Af(Y(s))ds) - —(f(Y(Nt)) - f(Y(0))). 



The first term on the right is a martingale and the second should go to zero, so if the 
martingale central limit theorem applies to the first, then Z N should converge. 

This paper addresses situations of the first type (V N =>■ Vo for a deterministic Vb, and we 
want to verify convergence of U N = r^iV^ — Vo)) in which both approaches are required. 
Specifically, the function F N giving the martingale, M , depends not only on Vq but also 
on another process (think V^it) = Vx(Nt)), so 

M N >\t) = V N (t) - Vf (0) - f t F N (V N (s),V 1 N (s))ds 

Jo 
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is a martingale, F N "averages" to F in the sense that 

-t 

?N { t rN { „\ t/JV/„\\ ipnrN/ 



(Vq*, Vf) is Markov with generator A N , and there exist H N such that A N H N w (F^ — F). 
(Note that if^v will be a vector of functions in the domain of A^, T>(Apj).) Assuming that 

ft 



M^(t) = (*)) - ^r(^r(0), ^ v (0)) - / A^^s), ^ v ( S ))ds 

Jo 

is a martingale, and again ignoring all the technicalities, we have 

r N (V N (t)-V (t)) (1.2) 
= nv(C(0) - V (0)) + r^M^f) - r N M N > 2 (t) 

+ J r N (F(V N (s)) - F(V (s)))ds 

+r N (H N (V N (t), V» (*)) - H N (V o N (0), Vf (0))) 

+r,v [\f n (V n (s), V»(a)) - F(V N (s)) - A N H N (V N (s), V»(s)))ds. 



If the last two terms on the right go to zero, the martingale terms converge 



r N M N ' 1 - r N M N ' 2 > M, 



and F is smooth, then we again should have U U satisfying (1.1). 

The work to be done to obtain theorems of this type is now clear. We need to identify 
F N and F, find an approximate solution to the Poisson equation A^Hn ~ F N — F, verify 
that the martingales satisfy the conditions of the martingale central limit theorem, and verify 
that the error terms (the last two terms in (1.2)) converge to zero. We will make this analysis 
more specific in stages. We are essentially considering situations in which the process 
is evolving on a faster time scale than Vq and "averages out" to give the convergence of 
Vq to Vq. But Vf* itself may evolve on more than one time scale. In the first stage of our 
development, we will replace Vy by (Vj , with and V^ 1 evolving on different (fast) 
time scales. Once the analysis for two fast time scales is carried out, the extension of the 
general results to more than two fast time scales should be clear. In the second stage, we 
consider multiply scaled, continuous-time Markov chains of a type that arise naturally in 
models of chemical reaction networks. For these models, many of the conditions simplify, 
but the notation becomes more complex. 



2 A central limit theorem for a system with determin- 
istic limit and three time scales 

We assume that V { N takes values in C E % i = 0,1, 2, and that E^ converges in the sense 
that there exists E, C M di such that E^ C Ej and for each compact K C R*, 

lim sup inf \x — y\ = 0. 
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We will refer to as the "generator" for the process V N = (V^, Vg 4 , V 2 N ), but all we 
require is that A N is a linear operator on some space V(A N ) of measurable functions on 
E N = Eg x Ef x Eg and that for h G V(A N ), 



h(V N (t))-h(V N (0))- / A N h(V"(s))ds 

Jo 

is a local martingale. 

We identify the time scales with two sequences of positive numbers and {r 2) Ar} 

and introduce a sequence of scaling parameters {tn} with the following properties. 

Condition 2.1 (Scaling parameters) The scaling parameters r^ — > oo and {rx,jv} and 
{r 2) Ar} are sequences of positive numbers satisfying 

limAr^no^- = 



rj,N 
T2,N 



lim^ooP = 0. (2 ' 1} 



Lo, Li, L2 will be linear operators defined on sufficiently large domains, T>(Lq) C M(Eo), 
X>(Li) C M(E x E x ), and £>(L 2 ) C M(E x E x x E 2 ), and taking values in M(E xEjX E 2 ). 
The requirements that determine what is meant by "sufficiently large" will become clear, 
but we will assume that the domains contain all C°° functions having compact support in 
the appropriate space. 



Condition 2.2 (Multiscale convergence) For each compact K C 

lim sup \A N h(v) - L h(v)\ =0, he T>(L Q ), 

lim sup I A N h(v) - Lih(v)\ = 0, h G 

N -*°°veKr\K N r l,N 

and 

lim sup I A N h(v) - L 2 h(v) \ = 0, he V(L 2 ). 

N ^°°v&Kr\E N r 2,N 



vdo+di+d,2 



Remark 2.3 Similar conditions are considered in Ethier and Nagylaki (1980). See also 



Ethier and Kurtz (1986), Section 1.7. 



There may be only two time-scales, in which case d 2 = 0, L 2 h = 0, and E = Eo x Ei (or 
equivalently, E 2 consists of a single point) in what follows. 

Condition 2.4 (Averaging condition) For each (vq,Vi) G E x E 1; there exists a unique 
Hv ,vi e V(E 2 ) such that J L 2 h(v ,Vi,v 2 )fi V0)Vl (dv 2 ) = for every h G V(L 2 )nB(E) . For each 
v G Eo, there exists a unique fj, Vo G P(Ei) such that f Lih(vo,Vi,v 2 )fi VQ)Vl (dv 2 )fi Vo (dvi) = 
for every h G V{L X ) fl ^(Eo x Ex). 

Remark 2.5 This condition ensures the uniqueness of the conditional equilibrium distribu- 
tion associated with the fast components. 
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With this condition in mind, we define 



Lih(v ,vx) = j L 1 h(vo,Vi,v 2 )fi vo , Vl (dv 2 ). 

For an K^-valued process Y, [Y] t will denote the matrix of covariations [Yj, Yj] t . 

Condition 2.6 (First convergence condition) There exist F N G M(E N , R d °) and F, G G 
C(E,M d °) such that 



M N,1 {t) = V N {t) - V N {0) - [ F N {V N {s))ds 

Jo 

is a local martingale, [V N ]t =>- 0, and for each compact K C E ; 



(2.2) 



lim sup \r N (F N (v) - F(v)) - G (v)\ = 0. (2.3) 

N ^°° v£KnE N 

Remark 2.7 This condition essentially implies L h = F ■ Vh for h G C^°(E ). 
Suppose that there exist hi G V(Li) d ° and h 2 , h 3 G V(L 2 ) d ° such that 

Li/ii(v ,vi) = J F(v , v 1 ,v 2 )fi V(hVl (dv 2 ) - J J F(vo,v 1 ,v 2 )fj IVOtVl (dv 2 )iJ, Vo (dv 1 ), 

L 2 h 2 (v ,vi,v 2 ) = F(v ,vi, v 2 ) - J F(v ,v 1 ,v 2 )ijLv ,v 1 (dv2), (2.4) 
L2h 3 (v ,vi,v 2 ) = Li/ii(u ,ui) - L 1 h 1 (v ,v 1 ,v 2 ). 
Define 

Fi(v ,vi)= / F(v ,vi,v 2 )fivo,vi{dv2), F(v ) = / / F(vo,v 1 ,v 2 )ii Vo>Vl (dv 2 )fiv (dvi), 



and 

H N = —h 1 + —(h 2 + h 3 ). (2.5) 

r l,N f2,N 

Note that for of this form 

A N H N « Lx/i! + L 2 (h 2 + h 3 ) = F -F, 



In what follows, if ^ does not have to be given by (2.5). That form simply suggests the 
possibility of finding H N with the desired properties. Specifically, we assume the existence 
of Hn G V(Ajsf) satisfying the following. 

Condition 2.8 (Second convergence condition) Assume that there exists G\ G C(E, R. d °) 
such that for each compact K C E, 

lim sup \r N (F(v)-F(v )-A N H N (v))-G 1 (v)\ = 0. (2.6) 

N ^°° v&KnE N 



Remark 2.9 The critical requirements for are (2.6), (2.10), and (2.11). In fact, because 

of the possibility of large fluctuations by and , even if h\, h 2 , and h 3 satisfying 

Condition 2.4 can be found, it may be necessary to define using a sequence of truncations 
of h\, hi, and h^. 



For Vq(0) G R d °, let V satisfy 



V (t) = V (0) + / F(V (s))ds 



(2.7) 



and define 



M N > 2 (t) = H N (V"(t))-H N (V N (0))- I A N H N (V"(s))ds. 



rJV 



rJV/ 



To identify the remaining conditions that are needed, we expand U N = r^{V^ 
U N (t) = U N (0)+r N (M N ' 1 (t)-M N ' 2 (t)) 

+r N [ (F(V N (s))-F{V (s)))ds 



N 



Vn 



(2- 



+r N / (F N (V N (s))-F(V N (s)))ds 
Jo 

+r N [ (F(V N (s)) - F(V N (s)) - A N H N (V N (s)))ds 
Jo 

+r N (H N (V N (t))-H N (V N (0))). 



The fourth term on the right is controlled by (2.3), the fifth by (2.6). (2.1) suggests that 
the sixth term goes to zero, but we will explicitly assume that. Assuming F is smooth, the 
third term is asymptotic to f* VF(Vo(s)) ■ U N (s)ds. 

That leaves the second term, and the following condition is needed for application of the 
martingale central limit theorem, Theorem A.l to this term. 

Condition 2.10 (Converegence of covariation) There exists G E C(E, M d ° xd °) such 
that for each t > 0, 

0, 



and 



lim E[swpr N \V "(s)-V "(s-) 
supr N H N (V N (s))^0, 

s<t 



(r N ) 2 [V N -H N oV%- f G(V N (s))ds 

Jo 



0. 



(2.9) 
(2.10) 

(2.11) 



Finally, we need a condition to ensure the relative compactness of the sequence. Let 
if) : E — > [l,oo) be locally bounded and satisfy lim^oo if)(v) = oo, and let denote the 
collection of continuous functions / satisfying 



sup 



< oo, lim sup 



0. 



For sequences of space-time random measures, the notion of convergence that we will use 
is that discussed in Kurtz (1992). 
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Lemma 2.11 Let V N be a sequence of K-valued processes, and define the occupation measure 



(2.12) 



(2.13) 



T N (Dx[0,t})= / l D (V N (s))ds. 
Jo 

Suppose that for each t > 

supE[ [ ?p{V N {s))ds} < oo. 

N Jo 

Then {T N } is relatvely compact, and ifY N ^-Y, then for fi, . . . f m E D^, 

( / h(V N (s))ds, f f m (V N (s))ds) ( / f\{v)Y{dv x [0, •]), . . . , / f m {v)T{dv x [0, •])) 
Jo Jo Je Je 

in C;rto[0, oo). 



Proof. Relative compactness of {T^} follows from Lemma 1.3 of Kurtz (1992). Relative 
compactness in Cr™ [0, oo) follows from relative compactness of each component. To see 
that for / e D^, the sequence X N = J f(V N (s))ds is relatively compact, it is enough 
to approximate the sequence by sequences known to be relatively compact. For e > 0, 
there exists a compact K e C E and C > 0, such that |/| < (Clx e + e)^. Define X^ = 
f Q 1 K e (V N ( s )) f (V N (s))ds . Note that X^ is Lipschitz with Lipschitz constant sup^g^ |/(v)|, 
so {X^} is relatively compact. For 5 > 0, 

supP{sup|X 7V (s) -Xf(s)| > 5} < UupE{ [ ^j{V N {s))ds], 
N s<t o n Jo 



and relative compactness of {X N } follows. (See Problem 3.11.18 of Ethier and Kurtz (1986).) 

Assuming that T N =>- T, the convergence of f Q f(V N (s))ds to J Ex , i f(v)T(dv x ds) 
follows by the same type of approximation. □ 



Condition 2.12 (Tightness) There exists a locally bounded ip : E — > [l,oo) satisfying 
lim^oo ip(v) = oo such that for each t > 0, 

supE[ [ tfj{V N {s))ds] < oo (2.14) 

N Jo 

and all of the following functions are in D^: sup^ \F N \, sup^ \r N (F N —F)\, sup N \r N (F—F— 
A N H N )\, \G\, su P7V \A N h\ for he V(L )nB(E ) , su PiV \^A N h\ fork e P(L 1 )n5(E xE 1 ) ; 

and sup N \^A N h\ for h e V(L 2 ) n B(E). 
Assuming the above conditions and defining 

G{vo) = J J G(v , v 1 ,v 2 )ijLvo,v 1 (dv2)nvo(dvi) (2.15) 

and similarly for Gq and G\, we have the following theorem. 
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Theorem 2.13 Under the above conditions, suppose that linijv^oo U N (0) = U(0), that F is 
continuously differentiate, and that the solution (necessarily unique) of (2. 7| ) exists for all 
time. Then for each t > 0, 

su P \V N (s)-V (s)\^0, 



s<t 



rjy(M ' — M ' ) =>■ M , where M has Gaussian, mean-zero, independent increments with 



E[M(t)M T (t)] = / G(V (s))ds 



(2.16) 



and U N =>- U satisfying 



U(t) = U(0) + M(t) + / (VF{V {s))U{s) + G {V {s))+G 1 {V {s)))ds. 



Assuming G = ao , we can write 
U(t) = U(0)+f a(V (s))dW(s)+ I (VF(V (s))U(s) + G (V (s)) + G 1 (V (s)))ds. (2.17) 



Remark 2.14 As noted above, the corresponding theorem for systems with two time-scales 
is obtained by assuming E 2 consists of a single point so L2/ = 0. 



Proof. Let be the occupation measure defined as in (2.12). Then by Lemma 2.11, {T^} 



is relatively compact. Assume, for simplicity that T. We will show that T is uniquely 

determined. 



Condition 2.6, (2.9), and the martingale central limit theorem, Theorem A.l, imply 
M N ' 1 0, and Lemma 2.11 then implies =>■ V£°, where 



V °°(t) = V (0) + / F(v)T(dv x ds). 

JEx [0,t] 



(2.18) 



Condition 2.12, the definition of L 2 , and Lemma 2.11 imply 
[ -(hi\ MO)) - / M.v//M M..s)W.s) = 

'Ex[0,t] 



2 N {h(V N (t)) - h{V N {0)) - / A N h{V N {s))ds) => I L 2 h{v)T{dv x ds) 



for every h E C^°(E). The uniform integrability implied by (2.14) implies that the limit is 



a continuous martingale with sample paths of finite variation and hence is identically zero. 



Condition 2.4 then implies (see Example 2.3 of Kurtz (1992)) that T can be written 



T(dv x ds) = fi V0 , Vl (dv2)T 0,1 (dv x dv\ x ds). 
A similar argument gives 

= / Lih(v)T(dv x ds) = / Lih(v Q , vi)T 0,1 (dv x dvi x ds), 

JEx[0,t] JE o xEix[0,t] 
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which implies 
t the convt 

Now (2.18b can be rewritten 



r 0,1 (<ifo x dv\ x ds) = fi VQ (dvi)T°(dvo x ds). 
But the convergence of V Q N to V^ 00 implies T°(dvo x ds) = 5v™( s )(dvo)ds. 



V~(t) = V (0)+ f F{V™{s))ds 
Jo 



and it follows that V °° 



Similarly, (2.11) now becomes 

(r N ) 2 [V N 



H N o V% 



G(V (s))ds, 



rN,2\ 





M as desired. 



(2.19) 



and it follows that r N (M N > 1 - M 

Finally, the uniform integrability implied by (2.14) and Condition 2.12 allows interchange 
of limits and integrals in the expansion of U N given in (2.8), and the convergence of U N to 
U follows. □ 



3 Diffusion approximation 



The functional central limit theorem, Theorem 



2.13 



suggests approximating V N by Vq + ^U. 
In turn, that observation and (2.17) suggest approximating V Q N by a diffusion process given 
by the Ito equation 



D N (t) = V N (0) + — / a(D N (s))dW(s 

rN Jo 
-t 



(3-1) 



/ ( F(D N (s)) + ±G (D N (s)) + ±-Gi(D»(s)) ) ds. 
o \ r N r N 



The approximation 



V N « 



D N = V + —U 

r N 



1 



N 



D 



N 



is, of course, justified by the Theorem 2.13 Justification for the approximation V c 
is less clear, since D N is not produced as a limit. Noting, however, that 

D N (t) = V N (0) + — f a(V (s))dW(s) 

rN Jo 

+ [ (f(V (s)) + —VF(V (s))U{s) + — G {V (s)) + — G^s)) ) ds, 
Jo V r N r N r N J 

assuming smoothness of F, G and G±, we see that r 2 N (D N — D N ) converges to U satisfying 



U(t) 



Va(V (s))U(s)dW(s) 



+ (vF(V (s))U(s) + ^U T (s)d 2 F(V (s))U(s) 



+{VG (V {s))+VG 1 (V (s)))U{s))ds 



and since the central limit theorem demonstrates that the fluctuations of V are of order 
0(r^- 1 ), we see that the difference between the two approximations D N and D N is negligible 
compared to these fluctuations. 



4 Markov chain models for chemical reactions 

A reaction network is a chemical system involving multiple reactions and chemical species. 
The kind of stochastic model for a network that we will consider treats the system as a 
continuous time Markov chain whose state X is a vector giving the number of molecules Xj 
of each species of type i G X present. Each reaction is modeled as a possible transition for 
the state. The model for the kth reaction, for each k 6 /C, is determined by a vector of 
inputs v k specifying the numbers of molecules of each chemical species that are consumed in 
the reaction, a vector of outputs v' k specifying the numbers of molecules of each species that 
are produced in the reaction, and a function of the state Xk{x) that gives the rate at which 
the reaction occurs as a function of the state. Specifically, if the kth reaction occurs at time 
t, the change in X is a vector of integer values (k = v' k — v k . 

Let Rk{t) denote the number of times that the kth reaction occurs by time t. Then R k is 
a counting process with intensity Afc(X(t)) (called the propensity in the chemical literature) 
and can be written as ^ 

Rk(t) = Y k ( [ X k (X(s))ds) } 
Jo 

where the Y k are independent unit Poisson processes. The state of the system at time t can 
be written as 

X(t) = X(0) +J2CkRk(t) = X(0) + Y^CkY k ( [ X k (X(s))ds) . 
k k Jo 

In the stochastic version of the law of mass action, the rate function is proportional to 
the number of ways of selecting the molecules that are consumed in the reaction, that is, 

A fe (x) = K k V ik \ W[ X }= K k W_Xi(Xi - l)---{Xi- V ik + 1). 

Of course, physically, \v k \ — Yli v ik is usually assumed to be less than or equal to two, but 
that does not play a significant role in the analysis that follows. 

A reaction network may exhibit behavior on multiple scales due to the fact that some 
species may be present in much greater abundance than others, and the rate functions may 



vary over several orders of magnitude. Following Kang and Kurtz (2012), we embed the 
model of interest in a sequence of models indexed by a scaling parameter N . The model 
of interest corresponds to a particular value of the scaling parameter N . For each species 
ieZ = {l,...,s}, we specify a parameter a« > and normalize the number of molecules by 
Nq* defining Z^°(t) = NQ 0i Xi(t) so that it is of 0(1). For each reaction k G /C, we specify 
another parameter (3 k and normalize the reaction rate constant as n' k = n k NQ k so that n k is 
of 0(1). One can observe this model on different time scales as well, by replacing t by tN^f, 
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for some 7 G R. The model then becomes a Markov chain on E^ = iV ai Z + x ■ ■ ■ x N as Z + 
which, when N = N Q , evolves according to 

Zfit) = Z?(0) + Y j N~ a X lk Y k ( f N^ a+h ^X^{Z N {s))ds) 
k Jo 

with 

If for some i, > and z^fc > 1, then Aj^ varies with iV but converges as iV — >■ 00. To simplify 
notation, we will write X k (z) rather than Aj^, but one should check that the iV-dependence 
is indeed negligible in the analysis that we do. Defining Ajv = diag(iV~ Ql , . . . , N~ aB ), so 
Z N = A N X, let 

A N f(z) = J2 N Pk X k (z)(f(z + A N ( k ) - f(z)), 

k 

where p k = v k ■ a + f3 k + 7. Since the change of time variable from t to tiV 7 is equivalent 
to scaling the generator by a factor of iV 7 , we initially take 7 to be zero. We subsequently 
consider the behaviour of Z N on different time-scales Z N (-iV 7 ). 

To be precise regarding the domain of An, note that because the jumps of Z N are 
uniformly bounded, if we define = inf{t : > r}, then for every continuous 

function /, 

f(Z N (t A r r N )) - f(Z N (0)) - / A N f(Z N (s))ds 

Jo 

is a martingale. 

For notational simplicity, assume that the aij satisfy < a± < • • • < a s , and let d Q > 
satisfy ctj = 0, i < d a , and a, > 0, i > d Q . 

To apply the results of Section[2j we identify r N , r ljA r, r 2)N from the reaction network and 
the parameters {f3 k } as follows. Let 

m 2 = max{p fe - a* : ( ik ^ 0}, 

and define r 2 ,N = N m2 . Then there exists a linear operator L 2 such that for each compact 
K c W, 

lim sup \—A N h(z) - L 2 h(z)\ =0, lie V{L 2 ) = C\R S ). 

Depending on the relationship between p k and ctj for 7^ and the time-scale parameter 7, 
the limiting operator L 2 is either the generator for a Markov chain, a differential operator, 
or a combination of the two, which would be the generator for a piecewise deterministic 
Markov process (PDMP). We classify the reactions by defining 

£2,0 = {k e K : p k = m 2 }, 

and 

^2,» = {k E K, : p k — oti = m 2 for some % with CKj > 0, Q k 7^ 0}. 
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For each k G K. 2j0 U 1C 2: , define 



( 2 , k = lim iV"^^ G Z. 

7V->oo 

Note that throughout the paper Q 2>k will denote the limiting reaction vector, not to be 
confused with the single matrix entry Q k . Then, for h G C 1 (IR S ) 

L 2 h(z)= *k(z){h(z + ( 2 , k )-h(z))+ J2 h(z)Vh(z) ■ ( 2>k . (4.1) 

Note that, although X k (z) depends on all species types, the dynamics defined by L 2 makes 
changes only due to reactions fC 2>0 U K 2 % . In other words, only the subnetwork defined by 
reactions K,ip U ¥i 2 ^ is relevant on the time-scale corresponding to 7 = —m 2 . If JC 2) » is 
empty, the process corresponding to L 2 is a Markov chain, and if K, 2 ,o is empty, the process 
is just the solution of an ordinary differential equation. If both are nonempty, the process is 



piecewise deterministic in the sense of Davis (1993) 



The process corresponding to L 2 can be obtained as the solution of 

V 2 (t) = V 2 (0) + V ( 2 , k Y k ( [ X k (V 2 (s))ds) + V C 2 , fc / X k (V 2 (s))ds, 

Jo Jo 



and assuming that V 2 does not hit infinity in finite time, Z N (-N~ m2 ) =>■ V 2 . 

The central limit theorem in Section [2] assumes that the state space is a product space 
and that the fast process "averages out" one component. The state space on which functions 



in the domain of L\ in Condition |2.2| are defined is such that every function on it is contained 
in the kernal of L 2 . In order to separate the state space in this way, we need to identify 
the combinations of species variables whose change on the fastest time-scale 7 = — m 2 is less 
than 0(1). This can be done with a change of basis of the original state space as follows. 

Let Sjc be a matrix whose columns are ( k , k G 7C for some subset 7C C /C. Then Sjc is 
the stoichiometric matrix associated with the reaction subnetwork K. For the species types 
whose behavior is discrete, Sk. gives the possible jumps, while for the species whose behavior 
evolves continuously, determines the possible paths. We will let TZ(Sjc) = spa.n{( k ,k G 
7C} C M s denote the range of S^, called the stoichiometric subspace of the chemical reaction 
subnetwork 7C, and we will let 



= {6eW: J^Cifc = OVA; G K} 



denote the null space of which is the othogonal complement of 1Z(Sk:)- For each initial 
value zo of the reaction system, z + 1Z(Sk.) defines the stoichiometric compatibility class 
of the system. Then both stochastically and deterministically evolving components of the 
system must remain in the stoichiometric compatibility class for all time t > 0. The linear 
combinations of the species 9 ■ X for 6 G Af^S^) are conserved quantities, that is, they are 
constant along the trajectories of the evolution of the reaction subnetwork 7C. 

On the time scale 7 = — m 2 , the fast subnetwork determined by L 2 has the stoichiometric 
matrix S 2 whose columns are {( 2)k , k G fC 2p U /C 2 ,.}- Define N(S 2 ) as above, and note that 
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9 ■ V-2, 9 G J\f(S%), are conserved quantities for the fast subnetwork, that is, 9 ■ V 2 (t) does not 
depend on t. Let s 2 denote the dimension of 1Z(S 2 ), and s[ = s — s 2 be the dimension of 
J\f(Sj). We now replace the natural state space of the process by J\f(Sj) x TZ(S 2 ), mapping 
the original processes onto this product space by the orthogonal projection ILy^T) x H-k(s 2 ), 
that is 

(v; N (t),v 2 N (t)) = (U M{s?) Z N (t),U n(S2) Z N (t)). 

Note that the original coordinates have different underlying state spaces N~ ai Z; however, 
the change of basis will combine only those coordinates with the same scaling parameter a*. 
To see that this is the case, note that by the definition of ( 2 ^ k , (2,ik 7^ and ( 2 j k 7^ 
implies a,i = Oij. It follows that there is a basis 9i,...,9 s > i for ftf(Sj) such that 9n ^ 
and 9ji ^ implies oti = aj, and we can take this basis to be orthonormal. We denote 
the common scaling parameter by ag r Let 6i be the matrix with rows 9j, . . . ,9j, so that 
(Qiz) T = (9 1 ■ z, . . . , 9 s ' i ■ z) T and the orthogonal projection is given by 



1=1 



On the next time scale we only need to consider the dynamics of the projection of 
the original process that is unaffected by the fast subnetwork V[ N = Hj^^ s t^K n X. Since 
n AT(5 2 T ) A JV = AjvIV^t), we have 

v[ N {t) = n mS T } z N (o) + a n n W) c*n(JV" f \ N k {z N {s))ds). 

Note that Hn(s 2 )(k is not necessarily equal to ( 2 ,k, nor is the other projection n^gr^ = 
(k — n^(5 2 )^ fc necessarily equal to ( k — ( 2jk . To identify the next time scale let 

mi = max{p fc - a 9l : 0, • ( k ^ 0} = max{p fc - on : (Il^^jCfc)* 7^ 0}, 

and define r^jy = N mi . Note that mi < m 2 . Then, there exists a linear operator L\ such 
that for each compact K C 

lim sup |— Ajv/i^) - Li/i(z)| = 0, (4.2) 

where h G P(Li) satisfies h{z) = f(0 1 -z,...,9 s[ - z) for / G C^R^). Define 

/Ci )0 = {k e K, : p k = mi, max |6>; • Cfe| > 0} 

and 

/Ci )# = {k e K, : p k — a dl = mi for some / with > 0, ^ • ( k ^ 0}. 

o — Ota . 

Let A® 1 = diag(iV^ a »i ,...,N < ), and for each fc G /Ci, U /Ci ># define 



d%= lim AT^-^A^Gia = lim N Pk ~ mi (N~ ae i9 1 AT ^ - Ckf 
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Then for h(z) = f{®iz) with / G C\W't) 

Uh{z)= \k(z)(f{e 1 z + tf ik )-f(e 1 z))+ h(z)vf(e lZ ).(l k . 

If V\ denotes the process corresponding to L\ then assuming that V\ does not hit infinity in 
finite time, V^ N (-N~ mi ) = U M(s t ) Z n (-N^ 1 ) V x . 

To separate the state space in terms of the next time scale (if there is one), define 

Ci lfc = Jim N^- m ^ N Yi msT) C k) 

in other words, £i,fc = QjCik * s embedded in the original space, and (f k — QiC k - On the 
time scale 7 = —mi, the subnetwork determined by L\ has the stoichiometric matrix S\ with 
columns {C,i, k k £ £1,0 U /Ci,»}. Define the subspace Af(Sf) as before, and let s\ denote the 
dimension of TZ(Si) and s = s[ — si be the dimension of Af(Sf). As before we need to map 
the processes onto this product space by the orthogonal projection IT^^T) x H^Si)- 
Since £ ljfc £ ■A/'(5'J) = span(6' 1 , . . . , S / ), we can assume that the ^ are selected so that 

ft(Si) = span(0 So+ i,...0 s >) = span(Ci,i,...,Ci, Sl )- 

Define 

so s i 

n = J2<WT = u ^st), n 1= = u n Sl) , and n 2 = u n{S2) . 

1=1 l=s +l 

On the next time scale we need only consider the projection HoZ N of the original process 
which is unaffected by either of the faster subnetworks. To identify the next time scale, let 

m = max{p fe - a Ql : 61 ■ ( k ^ 0, 1 < / < s } = max{p fe - on : (n Cfc)i ^ 0}, 

and define r ^ — N m °. Note that if 1 < s 2 , 1 < s±, 1 < s (s + «i + ^2 — s), m < mi < m 2 . 
Without loss of generality, we can assume that time is scaled so that m = 0. Then, there 
exists a linear operator L such that for each compact K C IR S °, 

lim sup \A N h{z) — L h(z)\ — 0, 

where /i G V(L ) satisfies = f(6 l -z,..., 9 S0 ■ z) for / G C\R S °). Define 

/C ,o = {k e K. : p k = 0,max |0* • Cfc| > 0} 

and 

/Co,, = {/c G /C : pit — ae, = for some / with > 0, 6\ ■ C, k 7^ 0, 1 < I < s }. 

As before, let O be the matrix with rows 6j,... 9j o , and let A® = diag(A^" Qe i , . . . , N~ at>3 « ), 
so that n = II^gT) = 0q ©o and for each k G /Co, U /Co,, define 

Cf'° = lim iV^A^0 o a = lim TV^^i ■ { k , . . . , A^o0 So • C,) T , 

iV— >oo TV — >oo 
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For h(z) = f(e z) with / G C^W ) 

L h(z)= M*)(/(e o s + tf°)-/(e *))+ E W/(Bo«) • Cf'°- 

To relate the above calculations to the results of Section |2l we assume that /Co o = so that 



L„M*) = £ A*(*)v/(e *) ■ #°. 

Let ^ = TZ W = (no^rixZ^ri^), so V" = (V N ^Vf ,V 2 N ) G Af(Sf) x ?e(5i) x 
TZ(S2), and note that T is invertible so that the intensities can be written as functions of v G 

Af(Sf)xn(S 1 )xn(S 2 ), that is A fc (T-^). Since U z = J2tx(0r z)B l and Yl lZ = Ezl 0+ i(#r 
the process Vq* = HoZ N is the embedding of QqZ n , and similarly (Vo V ,V 1 JV ) = 
(UqZ n , U±Z N ) is just the embedding of ®\Z N . Let Eo, Ei, and E2 denote the limit of 
the state spaces for , Vf* , and . 
The function F N in (2.2) is given by 



F N (v) = Y,N p *A%°\ k (T~ 1 v)e ( k 



(4.3) 



and 



F(v) = hm F N (v) = MT-Mcf' . 



To satisfy Condition 2.4 we will assume that L2 is such that for each (vq,vi) G Eo x Ei 
there exists a unique conditional equilibrium distribution yL Vom {dv 2 ) G V(E,2) for L2, Then 
Li/i(v ,«i) = / L 1 h( y vo,v 1 ,u 2 )fi vo , Vl (du2) is 



Li/i(>o,wi)= Ajfe(« ,ui)(/((«b,T;i) + Ci,fc) - /(uo,«i)) + 5^ A fc( w o,^i)V/(f ,fi) • Ci 



where A fc (f o,^i ) = / A fc (T l (t>o, fi, V2))Hvo ! v 1 (dv 2 ). For Condition 2.4 



(4.4) 

to be met, we also need 
to assume that for each v G E there exists a unique conditional equilibrium distribution 
^vo(dvi) G V(Ei) for L 1 . 

We further need to assume that there are functions hi G V(Li) : E x E x i— >■ Rl E °l and 
h,2,h,3 G T>{L<i) : E H- IRl E °l that solve the following Poisson equations: 



L 1 h l = F l -F, L 2 h 2 = F-F 1 , 



where 



Fi{v ,vx) 



L 2 h 3 = Lih x - L x h x , 
F{vQ,Vi,u 2 )^ VQ , Vl {du 2 ), F(v Q ) = / F 1 (v ,ui)fi vo (du 1 ) 



in order for Condition 2.8 to be met. We refer the reader to Glynn and Meyn (1996) for 



results on sufficient conditions for the existence of solutions to a Poisson equation for a 
general class of Markov processes. For the class of general piecewise deterministic processes 



see also Costa and Dufour (2003). For the examples considered in Section 5 we were able 
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to explicitly compute the desired functions. In general, however, explicit computation may 
not be possible, so results that ensure the existence of these functions may be useful. 

We now need to identify r^, which will be of the form = N p , for some < p < m x . 



Assuming that there is no cancellation among the terms in the sum in (4.3), for (2.3) to 
hold, we must have 



Then 



and 



p < max{a 9l - p k : 0, ■ ( k ^ 0, p k < a 8l , 1 < I < s }. 



6 t ■ G (v) = lim r N 0i ■ (F N (v) - F(v)) = V \ k {T^v)6 l ■ ( k 



(4.5) 



k:a Bl -p k =p 



so 



1=1 k:a Bl -p k =p 



Now let H N = r^hi + r 2 l N {h 2 + h%). To ensure that the limit in (2.6) exists, with 
reference to the definition of L 2 , we must have 



p < min{aj + m 2 - p k : ( ik 7^ °> a i + m 2 ~ Pk > 0} 



and 



p < min{2ai + m 2 - p k ■ Qk °, Qt f > 0}, 
and with reference to the definition of L\, we must have 



and 



p < mm{a ei + m x - p k : 0j ■ ( k ^ 0, a Bl + m 1 - p k > 0} 
p < min{2a e! + m x - p k : 9 t ■ ( k ^ 0, a 6l > 0}. 



(4.6) 
(4.7) 

(4.8) 
(4.9) 



Note that (4.5) implies the minimum in (4.8) and (4.9) only needs to be taken over so + 1 < 
/ < si. 

Assuming that hi, h 2 , and h s are sufficiently smooth, these assumptions insure that there 
exists Gi : E h-> Rl E °l 

/ A N A N \ 
Gi(v)= lim [r N ( L 2 ){h 2 + h 3 ) + r N ( L x )h x )= G l2 {v) + G n (v). 



r 2 M 



To identify G X2 , define 



,2,fc 



>1,N 



lim N p (N^- m *A N ( k - C 2 ,fc) 

AT— >oo 



and 



& >kij = lim N p+ »*- m *- ai - a i( ik ( kj . 

N—too 



{k G /C : 0/ ■ Cfc 7^ for some Z with = 0, m 2 — p k = p} 

{k e /C\/C 2) . : ■ (t ^ for some / with a e; > 0, m 2 — p k + = p} 
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Then setting h(z) = h^iTz) + h 3 (Tz), Guiv) = Hi 2 (T 1 v), where 

H 12 {z) = Y x k{z)Vh(z + C2,k)<2,k+ Y x k(z)S7h(z)-( 2tk 



+ Y ^k{z)\Y dz i dz N z )&M3 + Y X k(z)(h(z + ( 2 , k ) - h(z)) 



fee/cf 



Similarly, to identify Gn, define 



AT-j-oc 



iV— >oo 



and 



/Cf = {A; G /C : 61 ■ (k 7^ for some I with ae x = 0, mi — = p} 

^1 • = e : $2 " Cfe 7^ for some / with ae, > 0, mi — pk + ae t = p} 

Then G\\{v) = H n (T~ 1 v), where 

h u (z) = Y A fe (z)VM0i^ + Ci%)-C?, fc + Y ^vhie^-cU 

+ E wlY2 9 i d iM®iz)Zi,ku>+ Y h(z)(h 1 (e lZ + cl k )-h 1 (e 1 z)). 



We now need to identify G : E ->■ M^*™ satisfying ( |2.11[ ) in Condition 



2.10 



Let 



|Eo|x|E | 

^(t) = n(jv« f \ k {z N { s ))ds) 

Jo 

and ^(V*) = 60^^ - H N (V N ) = V N - H N {V N ). Then denoting z® 2 = zz T , 

N 2p [H N (V N )] t = Y N " P I (H N (V N (s-) + TA N ( k ) - H N (V N (s-)r 2 dR^(s), 
k Jo 

which is asymptotic to 

V N 2 ^ [\h n (V n (s-) + TA N ( k ) - H N (V N (s-)r 2 X k (Z N (s))ds. 
k Jo 

Taking the limit as iV -> 00 and integrating with respect to fi V0tVl (dv2) and fi Vo (dvi) then 
gives the value of G. 

5 Examples 

We now apply the central limit theorem to several examples of chemical reaction networks 
with multiple scales. 
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5.1 Three species viral model 



Ball, Kurtz, Popovic, and Rempala (2006) considered asymptotics for a model of an intra- 



cellular viral infection originally given in Srivastava, You, Summers, and Yin (2002) and 
studied further in Haseltine and Rawlingsj l 2002 ) . The model includes three time- varying 
species, the viral template, the viral genome, and the viral structural protein, involved in 
six reactions 



(1) 


T + stuff 




T + G 


(2) 


G 


m 


T 


(3) 


T + stuff 




T + S 


(4) 


T 


HA 





(5) 


S 







(6) 


G + S 




V 



whose reaction rates (propensities) are of mass-action kinetics form Xk{x) = n' k Yli x i ki with 
constants 



1 


1 


1 


K 2 


0.025 


2.5N 2/3 


K 3 


1000 


N 




0.25 


.25 




2 


2 


K 6 


7.5 x 10~ 6 


.75iV~ 5/3 . 



here expressed in terms of iVo = 1000. 

We denote T, G, S as species 1, 2, and 3, respectively, and let Xi(t) denote the number 
of molecules of species i in the system at time t. The stochastic model is 

Xi(t) = X 1 (0) + Y 2 ([ 0.025X 2 (s)ds) - F 4 ( f O^X^ds) 

Jo Jo 

X 2 (t) = X 2 (0)+Y 1 {[ X l {s)ds)-Y 2 {! 0.025X 2 (s)ds)-Y 6 ( [ 7.5 • lQ- 6 X 2 (s)X 3 (s)ds) 

Jo Jo Jo 

X 3 (t) = X 3 (0) + Y 3 ([ 1000X 1 (s)ds)-Y 5 (f 2X 3 (s)ds)-Y e ([ 7.5 • lO' 6 X 2 (s)X 3 (s)ds). 

Jo Jo Jo 

We take 

q;i = 0, a 2 = 2/3, 0:3 = 1. 
The scaling of the rate constants gives 



k 




0k 


Pk 


1 


1 








2 


2.5 


-2/3 





3 


1 


1 


1 


4 


.25 








5 


2 





1 


6 


.75 


-5/3 






Changing time t — > N 2 ^ 3 t, the normalized system becomes 

Z?(t) = Z? (0) + Y 2 { f N 2/3 2.5Z?{s)ds) - Y A { [ N 2/3 0.25Z? \s)ds) 

Jo Jo 

Z?(t) = Z?(0) + N-WY^ f N 2 ' 3 Z?(s)ds) - N~ 2 / 3 Y 2 ( f N 2 '*2.hZ» \s)ds) 

Jo Jo 

-N- 2 ' 3 Y,{ f N 2 / 3 0.75Z?(s)Z?(s)ds) 
Jo 

Z?(t) = Z?(0) + AT- 1 y 3 ( f N 5/3 Z^(s)ds)-N- 1 Y 5 ( [ N 5/3 2Z^ (s)ds) 

Jo Jo 

-N'^i [ N 2 / 3 0.75Z»(s)Z£(s)ds). 



We assume that the initial value for Z2 is chosen to satisfy Z 2 (0) = lim Z^iO) G (0, 00). 

iV-»oo 

In this model, there are only two time-scales, so we set 

, , , r 2 2 2 5 2,2 

mi = max{p fe - on : Qk T 0} = max{- - 0, 1, - - 1} = -, 

and we have r hN = N 2 / 3 . We have Ci,i = 0,Ci,2 = ei,Ci,3 = e 3 ,Ci,4 = -ei,Ci,5 = -e 3 ,Ci,6 = 
0. The operator L\ = lim7v-s>oo N~~ 2 / 3 A N is given by 

L x h{z) = \ 2 {z) (h(z + e x ) - h{z)) + A 4 (z) (h(z - e x ) - h(z)) + (A 3 (z) - \ B (z))d Z3 h(z) 

and note that for smooth h, 

N- 2 ' 3 A N h = L^ + O^N- 2 ' 3 ). (5.1) 

Functions h G ker(Li) are functions of the coordinate z 2 only, Ej = IZ(Si) = spanjei, e 3 }, 
and E = Af(Sf) = span{e 2 }. Taking h G V(L ) = C^Eq), 

L h(z) = lim A N h(z) = (Ai(z) - A 2 (z) - X 6 (z))d Z2 h(z 2 ). 

Setting V N = Z^ 1 and Vf = (Z± , Z^), the compensator for is 

F N (z) = \ 1 (z)-\ 2 (z)-\ 6 (z), 



so F(z) = F N (z) and G (s) = in Condition 2.6 



The process corresponding to L\ is piecewise deterministic with Z\ discrete and Z 3 contin- 



uous. For fixed z 2 , with reference to Condition 2.4 the conditional equilibrium distribution 
satisfies 



y 2.5z 2 (g{zi + 1, z 3 ) - g(zi, z 3 )) + ,25zi (#(zi - 1, z 3 ) - #(zi, z 3 )) 



(5.2) 



9o 

+(zi - 2z 3 ) — (zi,z 3 ) (j, Z2 (dzi,dz 3 ) = 0. 
(7Z3 
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Note that the marginal for Z\ is Poisson(102:2), so 

ziHz 2 (dzi, dz 3 ) = 10z 2 . 



Taking g(zi,z 3 ) = z 3 in (5.2), we see 

These calculations imply that the averaged value for the drift F is 

F(z 2 ) = / (Ai(V) - X 2 (z) - \ 6 (z))n Z2 (dzx,dz 3 ) = 7.5z 2 - 3.754, 



with WF(z 2 ) = 7.5 — 7.5^2. For the current example, we will see that F and G in (2.15) can 
be obtained without explicitly computing with fi Z2 . 



With reference to (2.4), we look for a solution h\ to the Poisson equation 
Lxh^z) = {z! - 2.5z 2 - 0.75z 2 z 3 ) - (7.5z 2 - 3.75^). 
Trying hi of the form hi(z) = ZxU\[z 2 ) + z 3 u 3 (z 2 ), we have 

L\h\{z) = u 1 (z 2 )(2.5z 2 - 0.25zi) + u 3 (z 2 )(zi - 2z 3 ) 



(5.3) 



and equating the factors multiplying z\ and z 3 , we get Ui(z 2 ) = 1.5z 2 — 4 and u 3 (z 2 ) = 0.375^2- 
Thus ht{z) = zi(1.5z2 - 4) + z 3 (0.375z 2 ) and H N {z) = N'^h^z). 



Since the solution of (5.3) is exact and (as we shall see) = N 1 ^ 3 , by (5.1), we have 
G\ = in Condition 2.8 With reference to Condition 2.10 ( |2.9 ) and (2.10) are immediate. 

The only restriction that remains to determine is the asymptotic behavior of the 
quadratic variation of Z o — H N (Z N ) = Z^ — N~ 2 ^ 3 hi(Z N ). Direct calculation shows that 
to get a nontrivial G in (2.11) we must take = N 1 ^ 3 . We then have 



nV 3 [z?-h n (z 



= J2 N ~ 2/3 (C2k + h l (Z N (s-))-h 1 (Z N (s-) + A N C k )) 2 dR^(s) 
k=i Jo 

« I Z?(s)ds+ [ (-1- 1.5Z?(s) + A) 2 2.5Z 2 N (s)ds 
Jo Jo 

+ f (l.5Z 2 N (s)-4) 2 .25Z?(s)ds + [ .75Z 2 N (s)Z 3 N (s)ds, 
Jo Jo 

where we observe that jumps by R 3 and R§ do not contribute to the limit. Dividing the 
equation for Z^ by iV 2//3 , we observe that 



Z?(s)ds 



10Z?(s)ds. 



Similarly, dividing the equation for Z% by N 2 ^ 3 we see that 



Z 3 (s)ds « - / Z?(s)ds n 



5Z 2 N (s)ds, 
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which in turn implies 

ft rt 



[ Z?(s)Z£(s)ds& [ hZ%{s) 2 ds. 
Jo Jo 



It follows that G(z 2 ) is 



10z 2 + (3 - l.hz 2 ) 2 2.hz 2 + (4 - l.hz 2 ) 2 2.hz 2 + 3.75^ 

= 72.5z 2 - 48.75*1 + 11.254- 



Let Z 2 be the solution of 



Z 2 (t) = Z 2 (0) + f (7.5Z 2 {s) - 3.75Z 2 {s))ds 
Jo 



and U N = N l ' 3 (Z* - Z 9 ). Then 



sup \Z${s) - Z 2 (s)\ =^ and U N => U 

s<t 



f ^72.hZ 2 (s) - 48.75Z 2 (s) 2 + 11.25Z 2 (s) W(s) + [ (7.5-7.5Z 2 {s))U{s)ds 
o Jo 



where, for W a standard Brownian motion, U satisfies 

u(t) = u(oy 

The corresponding diffusion approximation is 

D N (t) = (0) + N- 1 / 3 [ y/72.5D N (s) - A8.75D N (s) 2 + ll.25D N (s) 3 dW(s) 

Jo 

+ / (7.5D N (s) -3.75D N (s) 2 )ds 



We compare simulations for the original value of the amount of genome X 2 (-) with the 
approximations given by: the Gaussian approximation N 2 l 3 Z 2 (-N~ 2 l 3 ) J rN 1 l 3 U(-N~ 2 / 3 ), and 
the diffusion approximation N 2 I 3 D N (-N~ 2 I 3 ). For comparison we also give the deterministic 
value given by N 2 ^ 3 Z 2 (-N~ 2 /' 3 ). We use N = 1000 and a time interval on the scale 7 = 2/3. 
The initial values are set to -Xi(O) = X 3 (0) = 0,X 2 (0) = 10, and 500 realizations are 
performed for each of the three stochastic processes. Figure [T] shows the mean and one 
standard deviation above and below the mean for each of the three processes, and Figure [2] 
shows five trajectories for the three processes. 

For the diffusion process, these plots use only sample paths that hit one (= lOO/iv"^ 3 )) 
before they hit zero. For small initial values, the diffusion approximation does not give a 
good approximation of the probability of hitting zero (and hence absorbing at zero), before 
(for example) hitting one. Let 

= inf{t > : Z?(t) = or Z?(t) > 1] 

and 

= inf{* > : D N (t) = or D N (t) > 1}. 
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Figure 1: Mean and standard deviation of the amount of genome in the three species model 
(500 simulations with parameters JV = 1000, 7 = 2/3, Xi(0) = 0, X 2 (0) = 10, X 3 (0) = 0). 



It is shown in Ball et al. ([2006 ) that 



lim P{Z n (t") = 0\Z N (0) = N~ 2 ' z k} = A~ k 

iV— >00 

while a standard calculation for the diffusion process gives 



lim P{D N (rg) = 0\D N (0) = N^k} = e 



At 
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5.2 Michaelis-Menten enzyme model 

A basic model for an enzymatic reaction includes three time-varying species, the substrate, 
the free enzyme, and the substrate-bound enzyme, involved in three reactions 



(1) 
(2) 
(3) 



S + E ^ SE 
SE ^ S + E 
SE ^ P + E 



with mass-action kinetics and with rate constants such that k' 2 , k' 3 » k[. To be precise, let 



k' 3 = K3N, and k[ = K\ 



We denote E, S, P as species 1, 2, and 3, respectively, and let Xi(t) be the number of 
molecules of species i in the system at time t. Note that the total number of unbound and 
substrate-bound enzyme molecules is conserved, and we let M denote this amount. The 
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Figure 2: Five trajectories of the amount of genome in the three species model (same pa- 
rameters as in Figure [T]) . 
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stochastic model is 

X x {t) = X l {Q)-Y 1 {f K , 1 X 1 {s)X 2 {s)ds)+Y 2 {f k' 2 {M -Xi{s))ds) + Y 3 {f 4(M - Xi(s))ds) 

Jo Jo Jo 

X 2 {t) = X 2 {0)-Y 1 {[ K , 1 X l {s)X 2 {s)ds) +Y 2 { f k' 2 {M - Xi{s))ds) 

Jo Jo 

X 3 (t) = X 3 (0) + Y 3 ([ k' 3 (M -X 1 (s))ds). 

Jo 

If the initial amount of substrate is O(N) » M, then the normalizations of the species 
abundances are given by 

azi = 0, a 2 = 1, a 3 = 1 
and the scaling exponents for the rate constants are 

(3i = 0, f3 2 = 1, /3 3 = 1. 

The normalized system becomes 

Z?(t) = Z?(p)-Y 1 (f t NKiZ^(s)Z^(s)ds) + Y 2 ([ t NK 2 (M-Z?(s))ds) 

Jo Jo 

+Y 3 (f Nk 3 (M -Z N (s))ds) 
Jo 

Z*(t) = Z» r (0) - N- x Yi( [ NKiZ^(s)Z 2 N (s)ds) + N' 1 Y 2 ( [ Nk 2 (M - Z? \s))ds) 

Jo Jo 

Z*(t) = Zfty+N^Yaif ' NK 3 {M-Z?{s))ds). 

Jo 

Again, there are only two time scales with the fast time-scale mi = 1 giving r^jv = N. 
Then = — &i, Ci,2 = Ci,3 = e i; an d the operator L x is given by 

Lih(z) = Kiziz 2 {h(z - ei) - h{z)) + (« 2 + k 3 )(M - zi)(h(z + e x ) - h{z)), 

and for smooth h, 

N~ 1 A N h = Lih + 0(N- 1 ). (5.4) 

Functions h G ker(Li) are functions of coordinates z 2 and z 3 only. Thus Ei = {ziei : z\ = 
0, . . . , M} C K(Si) and E = A^-Sf ) = {(z 2 e 2 , z 3 e 3 ) : z 2 , z 3 > 0}. For h G £>(L ) = ^(Eq), 

LoM 2 ) = (k 2 (M - zi) - Kiz 1 z 2 )d Z2 h(z) + k 3 (M - z x )d Z3 h(z). 

Taking = {Z 2 , Z 3 ), the compensator for in (2.2) is 

F N (z) = (K 2 (M - Zi) - KiZiZ 2 , K 3 {M - Zl )f, 

so F(z) = F N (z) and G (z) = 0. 

On the fast time-scale, the process whose generator is Li is a Markov chain on Ei de- 
scribing the dynamics of an urn scheme with a total of M molecules, and for a fixed value of 
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z 2 , Z3, with transition rates K\z 2 for outflow and K2 + K3 for inflow. Its stationary distribution 
^z 2 ,z 3 ( z i) is binomial(M, p(z 2 )) for 

/ x ^2 + ^3 

= , : • 

So J (dzi) = Mp{z 2 ). 

This observation implies that the averaged value for the drift F is 

F(Z2, z 3 ) = (-M /^f , M -^g f = - K3 M(1 - P (z 2 )) ( \ 

with 

VP _ M «l«3(«2 + « 3 ) / 1 

(k 2 + % + k^) 2 V -1 
and we need to solve the Poisson equation 

LMz) = {k 2 {M - Zl ) - n lZl z 2 + M , k 3 (M - Zl ) - M ) T 

^2 + ^3 + K l z 2 K 2 + K 3 + K X Z 2 

for h x . Trying h\ of the form hi(z) = {z\U\{z 2 ), ZiU 2 (z 2 )) T , we have 

Lihi(z) = (~k 1 z 1 z 2 u 1 (z 2 ) + (k 2 +k 3 )(M-z 1 )u 1 (z 2 ), -hiz^u^z^ + ^+k^^M-z^u^)) 7 , 

and equating terms with the same power of z\, we get 111(22) — {^\Z 2 + k 2 )/(kiz 2 + k 2 + k 3 ) 
and u 2 (z 2 ) = k 3 /{k\z 2 + k 2 + k 3 ). Note that U\{z 2 ) + ^2(^2) = 1- Thus 

hl{z) = { (^+^+! 3 y {m + Z + k,/ = ZlMz ^ 1 - ni(Z2)f ' 

and H N (z) = N^h^z). 

Examining the quadratic variation of V N — H N o V N , we see that must be iV 1//2 , and 
by Q, it follows that G x = in Q. 

Finally, letting z® 2 = zz T ', 

N[V N -H N oV N ] t = N- 1 ^ [\Qo(k + h 1 (Z N (s-))-h 1 (Z N (s-) + A N ( k )y 2 dR%(s) 

k=i ^° 

T\ , ( U X {Z${ S )) \\«* KiZ N {s)z N {s)ds 



Oj \l- Ul (Z»(s)) 



<"1\ ( Ul (Z?(s)) ^ 2 



(<( (1 - ^l(^))) 2 -(1 - U^Z^s))? \ K 7 N, ]7 N ( ,, 

Jo V -(i - Mz 2 N (sW (i - Mz 2 N (s))f ) ^ ^ {s)ds 







-(l- Ml (Zf(s))) 2 (l-u^a)))' 
"/ ( n. (7N( B \\2 ,,)^N([.yl ) k 3 (M - Z^ (s))ds 



v -«i(^(*)) 2 «i(£?(*)) s 
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and averaging gives 

\im N[V N -H N oV N ] t = [ t G(Z(s))ds= t ( ^f* ( f\, ) *>, 

^oo L u U Jo Jo V -9( z 2{s)) g{Z 2 {s)) J 



where Z = (Z 2 , Z 3 ) satisfies 



Z(t) = Z(0) + f M «i«3^( a ) / -1\ ds 



o 



^2 + ^3 + «1^2(s) V 1 



and 



g{z 2 ) = M(l - «i(^ 2 )) 2 (/«ip(^2)2;2 + k 2 (1 - p{z 2 ))) 

+M M1 (Z 2 ) 2 /€3(1 -P(Z 2 )) 

Let 17" = iV 1 / 2 ^" _ Zzj z w _ Za )T_ Then 

sup|(Z^ (s) -Z 2 (s),Z^(s) -Z 3 (s))| =^0, and U N U 

s<t 

where U = (U 2 , U 3 ) T satisfies 

u(t) = u {0)+ f(-')^izM)dw (3)+ f ""tm^ 'fi'V 

Jo V 1 / io ( K 2 + ^3 + Kl^ 2 (s)) 2 V 1 / 

for W a standard scalar Brownian motion. 

The corresponding diffusion approximation is 

We compare simulations for 500 realizations of the original model with 500 realizations 
of the Gaussian approximation N Z 2 (-) + N^ 2 U 2 {-), N Z 3 (-) + N^ 2 U 3 {-) and of the diffusion 
approximation A^ -D 2 Y °(-), NqD^ ^). For comparison we also give the deterministic value 
given by N Z 2 (-), N Z 3 (-). We use N = 100 and a time interval on the scale 7 = 0. The 
initial values are set to Xi(0) = X a (0) = 0, X 2 (0) = 50 and M = 5, k\ = 0.1, k' 2 = 500, and 
k' 3 = 100. Figure [3] shows the mean and one standard deviation above and below the mean 
for each of the three processes, and Figure [4] shows five trajectories for the three processes. 
In this example, both Gaussian and diffusion approximations give good approximations for 
the means and the standard deviations of the pair of processes X 2 (-), X 3 (-). 



5.3 Another enzyme model 

Another model for an enzymatic reaction includes an additional form for the enzyme which 
cannot bind to the substrate. There are now four species, substrate, active enzyme, enzyme- 
substrate complex, and inactive enzyme, involved in five reactions 
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Figure 3: Mean and standard deviation of the amount of substrate in the Michaelis-Menten 
model (500 simulations with parameters JV = 100, 7 = 0, X x (0) = X 3 (0) = 0, X 2 (0) = 50, 
M = 5, k' x = 0.1, k' 2 = 500, 4 = 100). 

(1) S + E ^ SE 

(2) SE S + E 

(3) SE ^ P + E 

(4) F ^ E 

(5) E ^ F 

with mass-action kinetics and rate constants such that k[ = 0(1), k 2 , k' 3 = O(N), k' 4 , k' 5 = 
0(N 2 ) so that k' x = ki, k' 2 = k 2 X , k' 3 = k 3 N, k' 4 = k a N 2 , k' 5 = k 5 N 2 . 

We denote E, S, F as species 1, 2, and 3, respectively, and let Xi(t) be the number of 
molecules of species i in the system at time t. The total number M of active, inactive and 
substrate-bound enzyme molecules is conserved. The stochastic model is 

X 1 (t) = X 1 (0)-Y 1 ([ K , 1 X 1 (s)X 2 (s)ds)+Y 2 ([ 4(M - X 1 (s) - X 3 (s))ds) 

Jo Jo 

+Y 3 (f K ' 3 (M-X 1 (s)-X 3 (s))ds) + Y 4 ([ K 4 X 3 (s)ds)-Y 5 ([ ^X 1 (s)ds) 
Jo Jo Jo 

X 2 (t) = X 2 (0)-Y 1 (f K , 1 X 1 (s)X 2 (s)ds) + Y 2 ([ k' 2 (M-X 1 (s)-X 3 

Jo Jo 

X 3 (t) = X 3 (0)-Y 4 (f K , 4 X 3 (s)ds) + Y 5 ([ 4X X 

Jo Jo 



. 3 (s))ds) 
-i(s)ds). 
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Figure 4: Five trajectories of the amount of substrate in the Michaelis-Menten model (pa- 
rameters as in Figure [3]) . 
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If the initial amount of substrate is O(N) » M, then the scaling exponents for the species 
abundances are 

CKi = 0, a 2 — 1, «3 = 
and the scaling exponents for the rate constants are 

& = 0, (3 2 = 1, (3 3 = 1, (3 4 = 2, fa = 2. 

The normalized system becomes 

Z?(t) = Z?(0) - Yi( T N^{s)Z^{s)ds) + Y 2 {f Nk 2 (M - Z™{s) - Z%{s))ds) 

Jo Jo 

+F 3 ( / Nk z (M - Z?(s) - Z?(s))ds) + Y 4 ( [ N 2 K 4 Z^(s)ds) - Y 5 { [ N 2 K 5 Z?{s)ds) 
Jo Jo Jo 

ZN(t) = Z?(0) - N'^i f NK l Z^(s)Z 2 N (s)ds) + N~ 1 Y 2 ( f Nk 2 (M - Z?(s) - Z£{s))ds) 

Jo Jo 

Z?(t) = Z*(0)-Y 4 ([ N 2 K 4 Z 3 N (s)ds) + Y 5 ( [ N 2 K 5 Z?(s)ds). 

Jo Jo 

The fastest time-scale happens for m 2 = 2 and r 2 ^ = N 2 , with £ 2 ,4 = t\ — 63,^2,5 = 
— ei + e 3 . The operator L 2 is 

L 2 h(z) = k 4 z 3 [h(z + e 1 — e 3 ) - h(z)) + k 5 zi (h(z — e 1 + e 3 ) - h(z)) , 

with ker(L 2 ) consisting of functions of coordinates z 2 and z\ + z 3 only. To simplify our 
calculations we make a change of variables to (v , vi, v 2 ) = (z 2 , z\ + z 3 , z 3 ), so 

V Q N (t) = V N (0) - N-'Y.i /" jV«x(Vf (a) - V 2 N (s))V N (s)ds) + N^Y 2 { f Nk 2 (M - Vf (s))ds) 

Jo Jo 

V^it) = Vf (0) - Y 1 { f Nk^V^s) - V 2 N (s))V N (s)ds) + Y 2 ( f Nk 2 (M - Vf 

Jo Jo 

+Y 3 ([ NKsiM-Vfisfids) 
Jo 

V 2 N (t) = V 2 N (0)-Y 4 ([ t N 2 K 4 V 2 N (s)ds)+Y 5 ([ t N 2 K 5 (V 1 N (s)-V 2 N (s))ds), 

Jo Jo 

and in this system of variables £2,4 = e 2 , ^2,5 = — e 2 with the operator L 2 

L 2 h(v) = n 4 v 2 [h(v - e 2 ) - h(z)) + k 5 (vi - v 2 )(h(v + e 2 ) - h(v)). 

Functions h(v) G ker(L 2 ) are now functions of coordinates Vq,V\ only. Thus E 2 = TZ{S 2 ) = 
span{e 2 } and E] x E = Af(Sj) = span{ei,e }. 

The next time-scale corresponds to mi = 1, r ljN = N and = (0,-1), Ci,2 = Ci,3 — 
(0,1). Also 

Lih(v) = Kiv (v 1 -v 2 )(h((v ,vi) + (0,-i))-h(vo,v 1 ))) 

+(k 2 + « 3 )(M - «i) {h((v , Vl ) + (0, 1)) - h(v , Vl )) 
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with ker(Li) consisting of functions of v only. Thus Ei = TZ(S\) = spanjei} and Ei = 
Af(Sf) = span{e }. 
Finally, L is 

L Q h(v) = -Kiv (vi - v 2 )d Vo h(v ) + k 2 {M - vi)d V() h(v ). 

The Markov chain with generator L 2 is ergodic with, for a given value of (t> ,t>i), a 
stationary distribution /i V0jVl (dv 2 ) such that 



Po(vo,vi) = J v 2 fi V0 , Vl (dv 2 ) 



K-4 + K 5 

Thus the operator Li is 

L 1 h(v) = k iVo — (h({v , «i)+(0, -l))-/i(u , Ul)))+(K2+« 3 )(M-U 1 ) (/i((u , «i)+(0, l))-/i(«o, 

The Markov chain with generator Li is also ergodic with, for a given value of v , a stationary 
distribution p Vo (dv\) such that 

/ \ /" / 7 \ M(/t 4 + K 5 )(K 2 + K 3 ) 

Pi(^o) = / vin V0 {dvi) - 



KiK^Vq + (/C 4 + K 5 )(/t 2 + K3) 

and 

f \ f (a \ (a \ Mk 5 (k 2 + k 3 ) 

P2{Vo) = / v 2 fi V(hVl {dv 2 )/i V0 {dv 1 ) = 



K1K4V0 + (/C 4 + K 5 )(k 2 + «3) ' 

The compensator for the process Vj* is 

F N (y) = k 2 (M - Vl ) - «iu (ui - U2), 
so = F N (v) and G = 0. Averaging F gives 



^1(^0,^1) = y F(vo,v 1 ,v 2 )Hv ,v 1 (dv 2 ) = k 2 (M-vi) - kiv (vi -p (^o,^i)), 

and 

F(u ) = y ^1(^0,^1)^0(^1) = K 2(M - pi(v )) - K 1 v (p 1 (v ) - p 2 (v )) 

K^Vq «4(«2 + «3) 

= K 2 M — — K\MVq 



K1K4V0 + («4 + K 5 )(k,2 + K3) ' K1K4V0 + (/C4 + K 5 )(k 2 + K3) 



so 



/C1/C4UO + (/C4 + «5)(«2 + ^3) ' 

MK\K Z K4{K4 + K 5 )(k 2 + « 3 ) 



VF(v ) 



(k^Vq + («4 + « 5 )(re 2 + «3)) 



2 



Setting 



«i(^o) = —7 ; w ; V' "2(^0 J 



/C1/C4V0 + («2 + «3)(«4 + «5) " K4 + K5' 
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Figure 5: Mean and standard deviations of the amount of substrate in the three time-scale 
enzyme model (500 simulations with JV = 100, M = 5, 7 = 0, Xx(0) = 0, X 2 (0) = 50, 

x 3 (o) = 0, «i = 0.5, k' 2 = 500, 4 = 100, «4 = 4 = 5000). 

and 

U3{V0> = : Ui{Vq) 



the solutions to the Poisson equations Li/ii = F\ — F, L 2 h 2 = F — Fx, L 2 h 3 = Lih x — L\h\ 
are given by functions 

hx(v) = vxux(v ), h 2 (v) = -v 2 u 2 (vo), h 3 (v) = -v 2 u 3 (v ), 

andH N = ±hx + ^(h 2 + h 3 ). 

We again have G = 0, and G± = 0. We take r N = N 1 / 2 and observe that N~ 2 (h 2 + h 3 ) 
makes a negligible contribution to the quadratic variation. Consequently, we have 

5 , 

N[V N -H N oV N ] t » X> -1 / (Ck2 + hx(V N (s-))-hx(V N (s-)+TA N C k )) 2 dR^(s) 

k=i ^° 



(-1 + uxiVfis^KxVfiVfis) - V 2 N (s))ds 







+ / (l-u 1 (V* l (8))yK 2 (M-V»(a))d3 







+ / ux(V N (s)) 2 K 3 (M-Vx N (s))ds. 
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Figure 6: Five trajectories for the amount of substrate in the three time-scale enzyme model 
(same parameters as in Figure (5l) . 
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Hence 



G(v) 



K 3 {K 4 + K 5 ) 



+ 



( 



KxK^Vq + (k 2 + K 3 )(k 4 + K 5 

KlK 4 V + K 2 (k 4 + K 5 ) - - 



(kiv (vi - V 2 ) + K 2 (M - Vi)) 

(k 3 (M-Vi)) 



KKxK^Vq + («2 + K 3 )(k 4 + « 5 ). 

and £r(i>o) is obtained from the above function by replacing vi, v 2 by pi(fo) and p 2 {vo) giving 

M KiK^KiV^K^Ki + K 5 ) 2 (2K 2 + « 3 ) + (^lM^o + K 2 (k 4 + «5)) 2 ) 



G(« ) 



If Vn is the solution of 



V (t) = V (0) 



[k^Vq + (« 2 + «3)(«4 + K5)) 3 



M^K^V^s) 



/ K 1 K 4 V r (s) + (k 4 + ft 5 )(«2 + K3) 
then U N = AT 1 / 2 (Vq V — Vo) =>• £/ where, for W a standard Brownian motion, U satisfies 

U (s)as. 



U(t) = U(0)+ \ G(V (s))dW s - , 

Jo V Jo (^Vo(s) + (k 4 + k 5 )(k 2 + K 3 )f 

The corresponding diffusion approximation is 

Mk 1 k 3 k 4 D n (s 



D N (t) = Z? (0) + AT- 1 / 2 / v /G( J D JV (s))rfH/(,' 

Jo 



KxK 4 D n (s) + (k 4 + k 5 )(k 2 + K 3 



Finally, we compare simulations for 500 realizations of the original model X 2 with 

1 /2 

500 realizations of the the Gaussian approximation NqVq{-) + N Q [/(•) and the diffusion 
approximation NoD 2 °(-). For comparison we also give the deterministic value given by 
NqVq(-). We use Nq = 100, a time interval on the scale 7 = 0, and initial values are set to 
-X"i(0) = X 3 (0) = 0,X 2 (0) = 50 as in the previous example. Here the additional parameters 



are set to M = 5, k[ = 0.5, k' 2 = 500, k' 3 = 100, and k' 4 



5000. Figure 5 shows the 



mean and one standard deviation above and below the mean for each of the three processes, 
and Figure [6] five trajectories for the three processes. Again, both Gaussian and diffusion 
approximations give a good approximation for the mean and the standard deviation from 
the mean of X 2 (-). 

A Appendix 

A.l Martingale central limit theorem 



Various versions of the martingale central limit have been given by McLeish (1974), Rootzen 



ing version is from Ethier and Kurtz (1986), Theorem 7.1.4. 



(1977, 1980), Ganssler and Hausler (1979), and Rebolledo (1980) among others. The follow 
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Theorem A.l Let {M n } be a sequence of Mr-valued martingales. Suppose 

lim E[sup \M n (s) - M n (s-)\] = (A.l) 

n ^°° s<t 

and 

[MiMl) t ^ Ci>j (t), 

for all t > 0, where C = ((qj)) is deterministic and continuous. Then M n =>- M , where M 
is Gaussian with independent increments and E[M(t)M(t) T ] = C(t). 

Remark A. 2 Note that C(t) — C(s) is nonnegative definite fort > s > 0. If C is absolutely 
continuous, then the derivative will also be nonnegative definite and will have a nonnegative 
definite square root. Suppose C(t) = cr(t) 2 where a is symmetric. Then M can be written as 

M(t) = [ a(s)dW{s) 
Jo 

where W is d-dimensional standard Brownian motion. 

References 

Karen Ball, Thomas G. Kurtz, Lea Popovic, and Greg Rempala. Asymptotic analysis of 
multiscale approximations to reaction networks. Ann. Appl. Probab., 16(4):1925-1961, 
2006. ISSN 1050-5164. 

Rabi N. Bhattacharya. On the functional central limit theorem and the law of the iterated 
logarithm for Markov processes. Z. Wahrsch. Verw. Gebiete, 60(2):185-201, 1982. ISSN 
0044-3719. 

Oswaldo L. V. Costa and Frangois Dufour. On the Poisson equation for piecewise- 
deterministic Markov processes. SI AM J. Control Optim., 42(3):985-1001 (electronic), 
2003. ISSN 0363-0129. doi: 10.1137/S0363012901393523. URL |http : //dx . doi . orgTIoT 
I1137/S0363012901393523I 

M. H. A. Davis. Markov models and optimization, volume 49 of Monographs on Statistics 
and Applied Probability. Chapman & Hall, London, 1993. ISBN 0-412-31410-X. 

S. N. Ethier and Thomas Nagylaki. Diffusion approximations of Markov chains with two 
time scales and applications to population genetics. Adv. in Appl. Probab., 12(l):14-49, 
1980. ISSN 0001-8678. 

Stewart N. Ethier and Thomas G. Kurtz. Markov processes. Wiley Series in Probability and 
Mathematical Statistics: Probability and Mathematical Statistics. John Wiley & Sons 
Inc., New York, 1986. ISBN 0-471-08186-8. Characterization and convergence. 

Peter Ganssler and Erich Hausler. Remarks on the functional central limit theorem for 
martingales. Z. Wahrsch. Verw. Gebiete, 50(3):237-243, 1979. ISSN 0044-3719. doi: 
10.1007/BF00534147. URL |http : //dx . doi . org/10 . 1007/BF00534147"1 



34 



Peter W. Glynn and Sean P. Meyn. A Liapounov bound for solutions of the Poisson equation. 
Ann. Probab., 24(2):916-931, 1996. ISSN 0091-1798. doi: 10.1214/aop/1039639370. URL 



http : //dx . doi . org/10 . 1214/aop/1 039639370. 



Eric L. Haseltine and James B. Rawlings. Approximate simulation of coupled fast and slow 
reactions for stochastic chemical kinetics. J. Chem. Phys., 117(15) :6959-6969, 2002. 

Hye-Won Kang and Thomas G. Kurtz. Separation of time-scales and model reduction for 
stochastic reaction networks. Ann. Appl. Probab., 2012. (to appear). 

Thomas G. Kurtz. Limit theorems for sequences of jump Markov processes approximating 
ordinary differential processes. J. Appl. Probability, 8:344-356, 1971. ISSN 0021-9002. 

Thomas G. Kurtz. Strong approximation theorems for density dependent Markov chains. 
Stochastic Processes Appl, 6(3):223-240, 1977/78. 

Thomas G. Kurtz. Averaging for martingale problems and stochastic approximation. In 
Applied stochastic analysis (New Brunswick, NJ, 1991), volume 177 of Lecture Notes in 
Control and Inform. Sci., pages 186-209. Springer, Berlin, 1992. 

D. L. McLeish. Dependent central limit theorems and invariance principles. Ann. Probability, 
2:620-628, 1974. 

Rolando Rebolledo. Central limit theorems for local martingales. Z. Wahrsch. Verw. Gebiete, 
51(3):269-286, 1980. ISSN 0044-3719. doi: 10.1007/BF00587353. URL |http://dx.doi.| 
|org/10 . 1007/BF005873531 

Holger Rootzen. On the functional central limit theorem for martingales. Z. Wahrschein- 
lichkeitstheorie und Verw. Gebiete, 38(3):199-210, 1977. 

Holger Rootzen. On the functional central limit theorem for martingales. II. Z. Wahrsch. 
Verw. Gebiete, 51(l):79-93, 1980. ISSN 0044-3719. doi: 10.1007/BF00533819. URL 
|http : //dx . doi . org/10 . 1007/BF00533819| 

R. Srivastava, L. You, J. Summers, and J. Yin. Stochastic vs. deterministic modeling of 
intracellular viral kinetics. J. Theoret. Biol, 218(3):309-321, 2002. ISSN 0022-5193. 

N. G. van Kampen. A power series expansion of the master equation. Can. J. Phys., 39(4): 
551-567, 1961. doi: 10.1139/p61-056. 



35 



